## Plot raw data
## These results are reported in Appendix Figure A1
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

data4=data3

data4$Heat_Rate=data4$heat/data4$k_MWH
data4=subset(data4, k_MWH>0)


##############################################################################################################################
##############################################################################################################################
## Plot by treatment/control groups
data_base1=subset(data4, T<99)

data_base1$Group="Previous Years"
data_base1$Group[which(data_base1$Year==2018)]="Furlough Years"
agg=aggregate(data_base1[c("SO2_MASS..tons.","NOX_MASS..tons.",'AOD')], by=list(data_base1$Group,data_base1$T), FUN=mean, na.rm=T)
names(agg)[1:2]=c("Group","Time")
names(agg)[3:4]=c("SO2 (tons)","NOx (tons)")


setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "Raw_Trend_SO2.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(agg, aes(x = Time, y = `SO2 (tons)`, group=Group, color=Group),fill = "white") +
  geom_line(size = 1) +
  scale_color_manual(values = c("black","grey"))+
  scale_x_continuous(breaks=c(0, 25, 50, 75, 100, 125),
                     labels=c("Oct. 21", "Nov. 15", "Dec. 10", "Jan. 4", "Jan. 29","Feb. 23"))+
  geom_vline(xintercept=69, linetype="dashed", color='black', alpha=0.6)+
  geom_vline(xintercept=98, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=77, y=23, label="EPA furlough", size=8, color='black', alpha=0.6)+
  theme(legend.position="bottom", 
        panel.background=element_blank(), 
        panel.border = element_rect(fill = NA, color = "black"),legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))
dev.off()


setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "Raw_Trend_NOx.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(agg, aes(x = Time, y = `NOx (tons)`, group=Group, color=Group),fill = "white") +
  geom_line(size = 1) +
  scale_color_manual(values = c("black","grey"))+
  scale_x_continuous(breaks=c(0, 25, 50, 75, 100, 125),
                     labels=c("Oct. 21", "Nov. 15", "Dec. 10", "Jan. 4", "Jan. 29","Feb. 23"))+
  geom_vline(xintercept=69, linetype="dashed", color='black', alpha=0.6)+
  geom_vline(xintercept=98, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=77, y=17, label="EPA furlough", size=8, color='black', alpha=0.6)+
  theme(legend.position="bottom", 
        panel.background=element_blank(), 
        panel.border = element_rect(fill = NA, color = "black"),legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))
dev.off()

setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "Raw_Trend_AOD.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(agg, aes(x = Time, y = AOD, group=Group, color=Group),fill = "white") +
  geom_line(size = 1) +
  scale_color_manual(values = c("black","grey"))+
  scale_x_continuous(breaks=c(0, 25, 50, 75, 100, 125),
                     labels=c("Oct. 21", "Nov. 15", "Dec. 10", "Jan. 4", "Jan. 29","Feb. 23"))+
  geom_vline(xintercept=69, linetype="dashed", color='black', alpha=0.6)+
  geom_vline(xintercept=98, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=77, y=0.33, label="EPA furlough", size=8, color='black', alpha=0.6)+
  theme(legend.position="bottom", 
        panel.background=element_blank(), 
        panel.border = element_rect(fill = NA, color = "black"),legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))
dev.off()



